import cv2 as cv
import numpy as np


# 内容：统计米粒，读入灰度图片，增加灰度边框
image = cv.imread('/Users/apple/Desktop/imdata/rice.png', 0)
image = cv.copyMakeBorder(image, 3, 3, 3, 3, cv.BORDER_CONSTANT, value=63)
cv.imshow('1 - Original image', image)


# 构造5x5卷积核，5次腐蚀、5次膨胀，得到背景图
kernel = np.ones((5, 5), np.uint8)
erosion = cv.erode(image, kernel, iterations=5)
dilation = cv.dilate(erosion, kernel, iterations=5)


# 原图减去背景，得到米粒形状
backImg = dilation
riceImg = image - backImg
cv.imshow('2 - rice image', riceImg)


# 局部大津算法后，用开运算去噪点
ret, thresh = cv.threshold(riceImg, 0, 255, cv.THRESH_BINARY + cv.THRESH_OTSU)
element = cv.getStructuringElement(cv.MORPH_CROSS, (3, 3))
thresh = cv.morphologyEx(thresh, cv.MORPH_OPEN, element, iterations=1)
cv.imshow('3 - Threshold image', thresh)


# 轮廓检测
# GasImg = cv.GaussianBlur(thresh, (3, 3), 0)
# edge = cv.Canny(thresh, 50, 150)
# cv.imshow('4 - Canny image', edge)


# 检测闭合曲线、绘制轮廓、image转RGB图
contours, hierarchy = cv.findContours(thresh, cv.RETR_TREE, cv.CHAIN_APPROX_SIMPLE)
cv.drawContours(thresh, contours, -1, (120, 0, 0), 2)
image = cv.cvtColor(image, cv.COLOR_GRAY2BGR)


# 遍历米粒总数、得出平均面积、长度、3sigma米粒数量
count = 0
count3sigma = 0
area = []
length = []


for cont in contours:

    # 统计米粒数量、面积
    count += 1
    area.append(cv.contourArea(cont))


    # 得到矩形坐标长宽信息、统计长度
    rect = cv.boundingRect(cont)
    if rect[2] > rect[3]:
        length.append(rect[2])
    else:
        length.append(rect[3])


    # 逐一打印米粒信息
    print("{}-blob：{}".format(count, area[count-1]), end="  ")
    print("\tw:{} \th:{}".format(rect[2], rect[3]))

    # 绘制矩形
    cv.rectangle(image, rect, (0, 255, 0), 1)

    # 在米粒左上角写上编号，防止编号越出图片边界
    y = 10 if rect[1] < 10 else rect[1]
    cv.putText(image, str(count), (rect[0], y), cv.FONT_HERSHEY_COMPLEX, 0.4, (0, 0, 218), 1)


# 输出结果：平均面积
print("\n\n米粒平均面积：{:.2f}\t面积方差：{:.2f}".format(np.mean(area), np.var(area)))
print("米粒平均长度：{:.2f}\t长度方差：{:.2f}".format(np.mean(length), np.var(length)))


'''
# 计算3sigma米粒数量
upperArea = np.mean(area) + 3*np.var(area)
lowerArea = np.mean(area) - 3*np.var(area)

for cont in contours:
    if cv.contourArea(cont) > lowerArea & cv.contourArea(cont) < upperArea:
        count3sigma += 1

print("米粒总数量：{}\t3sigma内数量：{}".format(count, count3sigma))
'''

cv.imshow('Final image', image)


cv.waitKey()
cv.destroyAllWindows()

